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Despite a computational effort that scales exponentially with the system width, the Density 
Matrix Renormahzation Group (DMRG) method is one of the most powerful numerical methods for 
studying two dimensional quantum lattice systems. Reviewing past applications of DMRG in 2D 
demonstrates its success in treating a wide variety of problems, although it remains underutilized 
in this setting. We present techniques for performing cutting edge 2D DMRG studies including 
methods for ensuring convergence, extrapolating finite-size data and extracting gaps and excited 
states. Finally, we compare the current performance of a recently developed tensor network method 
to 2D DMRG. 



I. INTRODUCTION 

The behavior of quantum many-body lattice sys- 
tems depends strongly on their dimensionality. Mean 
field theory — where the behavior at each site is deter- 
mined self-consistently by the average influence of all its 
neighbors — works best when there are many neighbors, 
that is, in higher dimensions. Three dimensional systems 
can often be well understood using mean-field or semi- 
classical approaches, with quantum fluctuations acting as 
a minor correction.^ ^ In one dimension, mean field the- 
ory usually fails, and exotic behavior driven by strong 
quantum fluctuations, such as spin-charge separation, is 
the norm. By using powerful analytical and numerical 
approaches developed over the last few decades, it is pos- 
sible to determine the properties of ID systems to high 
accuracy. 

Two dimensional systems often have substantial quan- 
tum fluctuations, causing mean-field approaches to fail. 
Yet geometrical constraints are much more relaxed than 
in ID, allowing many more phases to exist. In some ways, 
geometry specifically selects 2D as being the most inter- 
esting: for example, in 2D a pair of particles can cir- 
cle one another, unlike in ID; but different numbers of 
revolutions (and their signs) are topologically distinct, 
unlike in three or more dimensions. This allows for the 
possibility — in 2D alone — of elementary excitations of a 
system's ground state which are anyons, particles that 
are neither bosons nor fermions. Currently, two dimen- 
sional systems are at the heart of quantum condensed 
matter physics, with a large fraction of researchers fo- 
cusing their efforts on key families of two dimensional 
systems, including the high temperature superconduct- 
ing cuprates,!^ quantum Hall systems^ and frustrated 
magnets, especially those that could host spin liquids.'^^ 

Most of the analytical and numerical techniques that 
often work in one or three dimensions fail in 2D. For ex- 
ample, analytic techniques based on the Bethe ansatz (an 
exact wavefunction for certain ID systems) or on confor- 
mal field theory are specific to ID.^^ Another example is 
quantum Monte Carlo, one of the most powerful numeri- 
cal methods. Classical Monte Carlo methods work in 
any dimension, but the path-integral formulation needed 
for quantum Monte Carlo introduces "probabilities" that 




FIG. 1: Results of DMRG simulations on 2D systems with ei- 
ther frustration or mobile fermions. Upper panel: local values 
of {S) for a triangular Heisenberg model. The 120° antiferro- 
magnetic order is "pinned" by applied magnetic fields on left 
and right edges; the top and bottom edges are connected to 
make a cylindrical geometry. Lower panel: spin and hole den- 
sities for a t- J cluster which exhibits both stripes and pairing. 
Cylindrical BCs are used, and pinning with both magnetic 
and proximity effect pair fields was used on the left and right 
open ends. 



can be negative; this is known as the sign problem.l^E^ 
Quantum Monte Carlo is of limited usefulness for sys- 
tems with a sign problem. Although the presence of the 
sign problem does not correspond strictly to dimensional- 
ity, ID systems with nearest neighbor interactions do not 
have a a sign problem, while 2D or 3D systems with either 
mobile fermions or frustrating spin interactions do have 
a sign problem. "Solving the sign problem" is sometimes 
called one of the great challenges of condensed matter 



2 



physics. Narrowly, solving the sign problem means fix- 
ing quantum Monte Carlo methods in some way, but we 
can also take a broader view: can we find any numeri- 
cal method that can solve broad classes of frustrated or 
fermionic quantum systems in 2D or 3D? 

The density matrix renormalization group (DMRG) is 
a numerical approach designed for one dimension and 
which has become the most powerful known numerical 
method in subject of this article is the appli- 

cation of DMRG to 2D. It is more difficult to use DMRG 
in 2D, and the results are much less accurate than in ID. 
Nevertheless, the lack of alternative approaches (when 
there is a sign problem) makes DMRG one of the most 
powerful current methods for certain 2D systems. It is 
a key goal of this article to present the tricks and tech- 
niques that are essential for the efficient treatment of 2D 
systems. 

Unfortunately, DMRG has been underutilized in this 
context. One reason may be the idea that since DMRG 
has a computational effort scaling exponentially with the 
width of the system, it is not useful. However, this same 
argument would imply that exact diagonalization (solv- 
ing for the ground state eigenvector utilizing the complete 
Hilbert space of a finite cluster) would not be useful even 
in ID ! In fa ct, exact diagonalization is often useful even 
in 2D. EES The exponential scaling of DMRG in 2D de- 
pends only upon the width, not the total number of sites 
like in exact diagonalization, and one often finds a fairly 
modest coefficient governing the exponential. One finds 
that even in sign-problem-free models, DMRG can make 
predictions of similar quality to Quantum Monte Carlo.^^ 

Figure 1 shows two examples of results from recent 
2D DMRG calculations. The system widths shown are 
close to the state of the art in terms of maximum sys- 
tem sizes reachable; but using the techniques discussed 



below in Section III, it is often possible to extrapolate 
such DMRG results to the infinite 2D limit quite reli- 
ably. DMRG makes up for its size limitations even fur- 
ther by providing access to the entire many-body wave- 
function, making it possible to compute essentially any 
ground state observable. 

To understand how DMRG works, observe that while 
exact diagonalization describes the wavefunction in a 
complete, exponentially large basis, some of these co- 
efficients and basis functions are much more significant 
than others. The idea of making a systematic approx- 
imation by truncating this basis is quite old, but in its 
simplest form it is not very effective for strongly corre- 
lated systems. In DMRG, one first rotates the basis so 
as to make the truncation much more accurate — in the 
rotated basis, only a few of the states are needed to rep- 
resent the ground state, while the rest can be discarded. 
A rotation involving the whole Hilbert space would be 
very inefficient; instead, DMRG uses many rotations fo- 
cused on a few sites at a time, generating global rotations 
via a sweeping procedure through all the sites of the lat- 
tice. The result is a wavefunction written in a particular 
form, called a "matrix product state" or MPS. An al- 



ternative viewpoint for DMRG is to start with the MPS 
form as a variational ansatz, then optimize all its coef- 
ficients. From this viewpoint, DMRG is an extremely 
efficient method for optimizing the coefficients.^^ The 
basis used by DMRG is the optimal basis within a cer- 
tain framework;^^ in recent years it has been understood 
to be a straightforward consequence of the Schmidt de- 
composition of quantum information. 

Because DMRG is so efficient, it has dominated nu- 
merical research into strongly correlated ID systems. 
One of the first applications of DMRG was calculat- 
ing the excitation gap in the S = 1 Heisenberg chain 
to very high accurac y a t a time when its existence re- 
mained controversial.^ Since then, not only has the 
original DMRG algorithm improved but its fiexibility 
has allowed many extensions. For example, DMRG is 
very useful for studying time-dependent Hamiltonian^^ 
and finite-temperature systems. Following the real- 
ization that DMRG is actually based upon matrix prod- 
uct states,'^ the method has become highly infiuential 
within the quantum information community, leading to 
the development of a new breed of algorithms based on 
tensor network states !^^*^ ' These developments have in 
turn accelerated the power and fiexibility of the original 
DMRG method. For example, it is now possible to use 
DMRG to reliably determine the ground state of infinite 
ID systems. 

While the hope is that tensor network state approaches 
will eventually overtake DMRG in two dimensional stud- 
ies (see Section IV), in the meantime it remains one of 
the most powerful and well controlled methods for sim- 
ulating models with a sign problem. DMRG can even 
be useful for simulating sign-problem free models since it 
provides full access to the many-body wavefunction.!^ 

In the following section we review past applications of 
DMRG to 2D systems. In Section |III] — a discussion of 
techniques for 2D — we assume a working knowledge of 
DMRG for ID systems; for a detailed introduction see, 
for example, Schollwock,!^^^ Hallberg,!^ NoaclP^ and 
White.i^ 



II. APPLICATIONS OF TWO DIMENSIONAL 
DMRG 

DMRG has been successfully used to study a wide va- 
riety of 2D systems - a testament to the fiexibility of the 
method. We briefiy survey some prior studies below, fo- 
cusing especially on those which break new ground by 
introducing techniques or characterizing especially inter- 
esting systems. The hope is to convey an understanding 
of what is possible with current methods and to motivate 
the investigation of systems that are good candidates for 
a DMRG study. 

One of the first applications of DMRG was in the area 
of magnetism.!^ Magnetic models, such as the Heisen- 
berg model, are ideally suited for DMRG because they 
can exhibit interesting phases with only modest entan- 
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glement and rarely suffer from large differences in en- 
ergy scales. Among the first 2D magnetic models to be 
studied with DMRG was the frustrated S = 1/2 Heisen- 
berg antiferromagnet on the CAVO lattice. This study 
helped establish the usefulness of DMRG for 2D sys- 
tems and introduced the wavefunction acceleration tech- 
nique, improving performance by up to two orders of 
magnitude. Since then, 2D DMRG has been used to 
study the frustr ated nea rest-neighbor Heisenberg model 
on the triangulai'^^'^^'^^' and kagome lattices^ as well 
as models frustrated by further neighbor or multi-spin 
interactions j ^^tt43 | Xwo dimensional DMRG is powerful 
enough to be a good option for studying unfrustrated 
spin models too.SSSSl One recent study obtained an esti- 
mate for the magnetization of the square lattice Heisen- 
berg model competitive with the best published Quan- 
tum Monte Carlo results.^ 

An area where 2D DMRG has been especially fruitful 
is in the search for quantum spin liquids — magnetic sys- 
tems that break no symmetries down to T = 0. Many 
realistic, one-dimensional spin models have disordered 
ground states, but the search for realistic 2D spin liquids 
remains a challenge.^ To date, DMRG has been used to 
identify tw o-dim ensional spin li quid p hases stabiliz ed by 
anisotropic further neighbop^^'^ and multi-spirP^'^ 
interactions, yet definitive evidence of a short-range, 
isotropic spin model whose ground state breaks no sym- 
metries has been lacking. However, DMRG simulations 
of the Heisenberg antiferromagnet on the kagome lattice 
now show strong evidence for a spin liquid ground state.^ 
2D DMRG has also been used recently to study pertur- 
bations to the 'unrealistic' Kitaev honeycomb model that 
may actually describe a certain limit of the layered mag- 
net NasIrOs^^ 

Another set of systems studied extensively with 2D 
DMRG are the doped t - J and Hubbard models. DMRG 
is an attractive option for simulating models of this 
type because it can deal with arbitrary doping in an 
unbiased way. A common thread running through the 
DMRG literature on the t — J and Hubbard models 
is the formation of stripes - charge density waves sep- 
arating regions of phase-shifted antiferromagnetic or- 
der. After stripes were found to occur spontaneously 
in four- leg t — J ladders,^ they were also found to oc- 
cur in systems up to width eight .^^ Follow-up studies 
have examined properties of stripes more closely, such 
as their doping and interaction energ y/ and the effects 
of further-neighbor interactionJ^^^^ and anisotropic^^ 
likely present in real materials. It was suggested early 
on that the stripes observed experimentally are due to 
frustrated phase separation,*^ but evidence from DMRG 
favors the idea that stripes are caused by local compe- 
tition between hopping and exchange, without the need 
for long range Coulomb interactions.^^ Two-dimensional 
DMRG results may also shed light on the pairing mech- 
anism of the cuprates. 

One study which especially highlighted the reliability 
of DMRG results for 2D t — J models was Chernyshev 



et al.'s comparison of 2D DMRG results to predictions 
from a self-consistent Green function theory of holes in an 
antiferromagnetic background. Working with a modified 
t — J model where J is replaced by an Ising Jz interac- 
tion, the two methods show remarkable correspondence 
in their predictions of energies and hole d istributions, es- 
pecially for the case of a sing le hole."^^ Further DMRG 
investigations of the t — J model have examined the ex- 
istence of checkerboard order, edge states of holes in 
nanosystems^^ and the competition between stripes and 
pairing. DMRG has also been used to demonstrate a 
striped phase in the Hubbard model on systems up to 
width six.lSDEH 

Another class of systems well suited for 2D DMRG 
are frustrated bosonic models, although there have been 
few DMRG studies of these systems until very recently. 
Bosonic models can be simple to work with numerically, 
yet are expected to exhibit new kinds of phases and ex- 
otic orders not found in one dimension. They also have 
great potential for controlled experimental realizations 
with cold atoms. One topic in bosonic systems that has 
attracted great interest is the possibility of a supersolid, 
a phase with simultaneous charge density wave and su- 
perfluid order. A supersolid phase was shown to exist 
in the unfrustrated triangular lattice model, but the or- 
der turned out to be quite weak. However, Jiang et al. 
have demonstrated using 2D DMRG that adding further 
neighbor hopping to the model significantly stabilizes su- 
persolid order 

A rather different kind of bosonic phase predicted to 
exist in 2D is known as a bose metal. This phase is 
not captured by a conventional order parameter, but in- 
stead by a pattern of correlations associated with an en- 
tire surface of gapless modes (a 'bose surface'). Using 
DMRG, it has proved possible to identify remnants of the 
full 2D surface in quasi- ID models with ring-exchange 
interactions on two-leg^ and four-leg^^ ladders. Moti- 
vated by the original construction of the bose metal in 
terms of fermionic partons, DMRG calculations have also 
been used to show that a frustrated 2D fermionic model 
could enter a Cooper-pair bose metal phase. Recently, 
a rather different type of ordering was observed for a 
fermionic system on the kagome lattice, where 2D DMRG 
was used to study a metal-insulator transition.!^ 

One final study that does not fit neatly into the above 
categories but certainly bears mentioning is the investi- 
gation by Jeckelmann and White of the Holstein model, 
which describes a single fermion interacting with bosonic 
lattice vibrations or phonons.^^ An effective way of deal- 
ing with the large number of boson species in this model 
was to split each site into many smaller ones having fewer 
degrees of freedom. Using this technique it proved pos- 
sible to simulate cylindrical systems with widths up to 
twenty! 



III. TECHNIQUES FOR TWO DIMENSIONAL 
DMRG 

When working at the frontier of current numerical ca- 
pabihties, it is important to ensure that individual sim- 
ulations give trustworthy results and then combine these 
results correctly to build an accurate picture of a 2D 
model. For two dimensional DMRG, the main obstacle 
to overcome is that the number of states kept must be 
increased exponentially with the width of the system to 
maintain a constant accuracy/^ In practice, this puts an 
upper bound on the system sizes that can be simulated. 
In addition, one should use open or cylindrical bound- 
ary conditions, as opposed to fully periodic, in order to 
avoid squaring the number of states required for a given 
accuracy (see below). 

Although the first restriction is unfortunate, the use of 
cylindrical boundary conditions is hardly the drawback it 
is often portrayed to be. In fact, cylindrical boundaries 
provide a significant degree of control over simulations 
and can be very useful for inferring properties of the full 
2D system. The limitations of two dimensional DMRG 
are further mitigated by the enormous fiexibility of the 
DMRG approach. DMRG gives full access to the many- 
body wavefunction. This means, for instance, that the 
ground state computed for one Hamiltonian can be used 
as input for simulations with a different Hamiltonian, 
inviting a range of techniques such as biasing initial states 
to detect symmetry breaking or changing a Hamiltonian 
in mid-simulation to locate phase boundaries. 

In what follows we discuss these techniques and others 
that have been successfully used in cutting-edge calcula- 
tions. The hope is to equip practitioners to pursue new 
state-of-the-art 2D DMRG calculations and to motivate 
them to create their own techniques, extending the reach 
of current numerical methods. 



A. Ensuring Ground State Convergence 

A successful DMRG study of a two dimensional sys- 
tem involves multiple independent calculations for vari- 
ous system sizes over a range of parameters. In order to 
deduce the correct behavior of the 2D system from these 
results, one must guarantee that each ground state cal- 
culation is well understood and as accurate as possible. 
DMRG can fail to find the true ground state for two basic 
reasons. The first is that the number of states kept after 
each truncation of the wavefunction may be too small to 
represent the wavefunction accurately. The second issue 
is more subtle and arises because while DMRG uses exact 
diagonalization locally, it is globally a variational method 
and can therefore get stuck in a metastable state. 

There are a number of ways to ensure that one does 
not over-truncate the wavefunction and that a sufficient 
number of states m are being kept. The most impor- 
tant is to avoid using fully periodic boundary conditions, 
as they require states to represent the same wave- 
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FIG. 2: Interaction bonds for a nearest-neighbor model with 
cylindrical boundary conditions on a three leg ladder. The 
thicker bonds indicate interactions which remain nearest- 
neighbor along the ID path used by DMRG while thinner 
bonds are treated as further neighbor interactions. 



function requiring only m states with open or cylindrical 
boundary conditions. By cylindrical boundary condi- 
tions we mean open boundary conditions along the larger 
lattice direction (the 'length', or x direction) and peri- 
odic boundary conditions along the smaller direction (the 
'width'). A typical setup using cylindrical boundary con- 
ditions is shown in Figure [2j 

After a calculation is under way, the number of states 
kept should be increased systematically until the energy 
and any other observables of interest converge to a spec- 
ified tolerance. This can be done by directly adjusting 
m, or by fixing the maximum truncation error at a step 
and have the code determine m on the fiy. An optimal 
approach may involve taking both a maximum (and min- 
imum) m into account as well as a requested truncation 
error. A minimum m is important because the reported 
truncation error early on in the calculation could be very 
inaccurate, resulting in slow convergence because of too 
small an m. 

Metastability issues are more difficult to deal with sys- 
tematically. Practitioners must use a number of checks 
in addition to physical insight to verify that the wave- 
function has reached the true ground state. For many 
models, specific properties of the ground state are known 
in advance. Additionally, if one uses a DMRG code that 
represents wavefunctions as a matrix product state, then 
by representing the Hamiltonian as a matrix product op- 
erator it is possible to efficiently compute the variance 
o'H = {i^\H'^\i^) — E'^ and confirm that is an eigen- 
state. 

There is a general technique which helps DMRG avoid 
being stuck in metastable states. This technique is par- 
ticularly helpful in making sure the MPS wavefunction 
includes correlations coming from terms in the Hamil- 
tonian connecting distant sites in the MPS path. For 
example, a nearest neighbor hopping in a 2D strip may 
connect sites which are the transverse width apart in the 
MPS path. These correlations may have trouble getting 
started, since the extra states needed to allow the hop- 
ping at one of the two sites may not help lower the en- 
ergy unless the extra states at the other site are already 
present. A special "noise" term can be added to the 
density matrix at each step which takes into account all 
terms connecting the left and right blocks. This ex- 
tra noise helps such long range correlations get started; 
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generally the noise is turned off in later sweeps. This 
technique is particularly important if, for some reason, 
fully periodic boundary conditions must be used. 

A key way to avoid metastability is to begin with a 
wavefunction that is already close to the true ground 
state. Good initial wavefunctions may be hard to come 
by for an unfamiliar model; here careful calculations on 
smaller system sizes or at smaller values of m where 
DMRG has better control can be an excellent guide. Such 
calculations allow one to identify the dominant correla- 
tions within the ground state which may otherwise be 
obscured by strong fluctuations on larger lattices or for 
larger m. 

For a system which is expected to have a conventional 
symmetry breaking ground state, such as an antiferro- 
magnet on a bipartite lattice, a Neel state may be a 
sufficiently good starting point. For systems with more 
subtle order, the initial wavefunction can be produced by 
starting with a lower symmetry or 'pinned' Hamiltonian. 
For example, if a system is expected to have a ground 
state with valence bond solid order, the Hamiltonian can 
be modified by adding pinning fields A Si ■ Sj for each pair 
of sites i, j connected by a valence bond. Then, after a 
few sweeps A can be gradually tuned to zero allowing the 
system to relax to its true ground state. 

Using an initial state or a pinning field can also be 
helpful for ruling out hypothesized properties of a model. 
If DMRG restores a symmetry explicitly broken by the 
initial state, one has strong evidence against that par- 
ticular ordering scenario. This method has been used to 
rule out a type of checkerboard order for the t—J modeP^ 
and more recently as evidence against a particular type 
of valence bond solid order for the kagome Heisenberg 
antiferromagnet ^ 

When dealing with a complex lattice or a phase with 
a large unit cell, one way to deal with metastability is- 
sues and minimize the number of states needed is to ex- 
periment with multiple DMRG paths. For a fixed value 
of m, DMRG is better able to capture entanglement on 
Hamiltonian bonds that remain nearest-neighbor when 
mapped to ID. Choosing the DMRG path judiciously 
can even permit complex initial wavefunctions such as 
valence bond solids to be represented exactly with only 
a small value of m. Having the ability to reproduce 
the same ground state with different DMRG paths can 
also provide strong evidence that one has found the true 
ground state and not a metastable solution. 

After having gained a good understanding of smaller 
systems, one wants to push DMRG calculations up to the 
largest accessible widths. At these widths, there is less 
control, so in order to produce accurate results it is very 
useful to extrapolate from more controlled limits. For 
DMRG, a natural extrapolation parameter is the trunca- 
tion error £ (the sum of discarded density matrix eigen- 
values). The energy has long been extrapolated to zero 
truncation error, where normally a linear extrapolation 
of E versus £ is best. Remarkably, within a DMRG calcu- 
lation local measurements performed on the two central 




FIG. 3: Results of a DMRG calculation for the Heisenberg 
model on a 16x8 cylinder with antiferromagnetic order pinned 
at the open boundaries. To work in the strong pinning limit, 
it is useful to imagine the finite system embedded within a 
larger system acted on by an infinitely strong field (shown 
here as the shaded regions) . The pinning fields at the physical 
edges are determined by the Hamiltonian bonds connecting 
the real and fictitious system. 

sites at each step also have errors varying linearly with £\ 
This is one reason why measuring local quantities, per- 
haps in response to a perturbation, is usually preferred to 
correlation functions, whose error varies as e^/^.^^ Effi- 
cient extrapolations can be performed using results from 
a single DMRG calculation with increasing m, but it is 
important to repeat each m for two full sweeps (and to ex- 
trapolate using the last of the four half-sweeps) to ensure 
that the calculated £ is consistent enough for extrapola- 
tion. 

The flexibility of DMRG even allows the use of other 
extrapolation parameters that may work better than the 
truncation error in certain cases. For example, one can 
add a small perturbation A H' to the Hamiltonian and ex- 
trapolate the energy in A. For this approach to work well, 
the ground states of the perturbed Hamiltonian should 
be less entangled than the true ground state. Further- 
more, by choosing H' to have a vanishing expectation 
value with respect to the ground state of the first 
derivative of the energy with A can be tuned to zero, 
increasing the accuracy of the extrapolation.^ 

B. Working Around Finite Size Limitations 

A variety of approaches can be taken to predict bulk 
2D behavior from sets of finite systems - here we dis- 
cuss some that are particularly useful in the context 
of DMRG. Most of these approaches utilize cylindrical 
boundaries. On the two open edges of the cylinder, one 
is free to apply local fields ("pinning" it), or to perturb 
it in other ways, in order to make the bulk represent 
2D most accurately. A favorable side-effect of applying 
a boundary pinning field may be a reduction of entan- 
glement, improving the DMRG convergence. For exam- 
ple, an antiferromagnet on a finite system typically has 
a singlet ground state, but one may regard it as a su- 
perposition of antiferromagnetically ordered states with 
different directions for the order parameter. Pinning can 
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select one order parameter direction, reducing the com- 
plexity and entanglement of the state, while simultane- 
ously representing the broken-symmetry 2D properties 
more faithfully.^^ 

It can be very useful to choose strongly pinned edges. 
To choose this type of pinning, it is helpful to think of 
the finite cylindrical system as part of an infinite cylin- 
der. A fictitious infinitely strong pinning field is applied 
to the fictitious system outside the finite part. The ficti- 
tious part then acts on the boundary of the real system 
through the Hamiltonian terms connecting the two. Fig- 
ure |3] shows the results of a ground state calculation for 
the 5 = 1/2 Heisenberg model on a cylinder of width 8 
and length 16. Two rows of the fictitious system are also 
shown at each end of the cylinder. The magnetization 
is enhanced near the pinned sites, then approaches the 
2D bulk value at the center (if the aspect ratio is chosen 
properly as discussed below). 

After gaining a good understanding of the system's 
ground state properties, the next step is to simulate mul- 
tiple system sizes to build up an accurate picture of the 
thermodynamic limit. Because DMRG scales much more 
favorably with the length of a system than with its width, 
it is a good idea to group results at fixed width together. 
Naively then, it would seem that the best approach is to 
simulate the longest possible system at each width. But 
there are much more efficient ways to proceed. 

For a fixed width, the energy per site of a length 
system approaches the infinite system value with an er- 
ror proportional to l/L^ because of a constant term in 
the energy from the open ends. This slow convergence re- 
duces the efficiency and accuracy of direct extrapolations. 
Therefore, it is convenient to determine bulk cylinder en- 
ergies by subtracting the energies of different finite-length 
cylinders. As illustrated in Figure [4j subtracting cancels 
edge effects, leaving only the bulk energy of the larger 
system, which rapidly converges to the infinite value as 
a function of L^- The convergence rate depends on the 
bulk correlation length(s). If the correlation lengths are 
infinite, one can extrapolate the subtracted energies with 
a polynomial for the highest accuracy. 

As an example of the subtraction method. Figure |5] 
show results for the square lattice 5* = 1/2 Heisenberg 
model on 6 leg ladders of length = 8 through Lx = 20. 
The bare energies of the finite ladders have a strong 
l/Lx dependence, but the subtracted energies rapidly 
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FIG. 5: Subtraction method for estimating the energy per 
site of the 5* = 1/2 Heisenberg model on an infinite cylinder 
of width Ly = 6. The upper points (squares) are the energy 
per site of finite systems with Lx = 8, 10, 12, 16 and 20 - er- 
ror bars were smaller than the symbol size. The lower points 
(triangles) are found by subtracting the energy of the length 
Lx/2 system from the length Lx system and dividing by the 
number of extra sites of the larger system. The energy esti- 
mate Co = —0.67172(7) for the infinite cylinder is taken from 
the best subtracted data point. 
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FIG. 6: Extrapolation of energy estimates for infinite cylin- 
ders of the S = 1/2 Heisenberg model with J = 1. The 
energy densities for widths Ly = 6, 8, 10 and 12 are fit to a 
cubic polynomial a -\- b Ly^ (see Ref. 73 ), giving an energy 
estimate eo = —0.6694(3) for the 2D system. 
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FIG. 4: Properties at open edges are not typical of the bulk 
model and may include pinning fields or extra sites designed 
to accommodate a certain ground state. Subtracting energies 
calculated for two finite length system cancels the leading edge 
effects giving values that rapidly converge to the properties 
of the infinite cylinder. 



converge to a constant value and require no fine tuning of 
the boundary conditions. While one has to extrapolate 
the single-cylinder data, the subtracted energies converge 
quickly enough that the best point can be used as the in- 
finite cylinder energy estimate. 

Finally, after obtaining the bulk energies of infinite 
cylinders with various widths Ly, one can extrapolate 
in Ly to estimate the 2D energy. For the case of pe- 
riodic boundary conditions — applicable for the infinite 
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FIG. 7: Ground states of the S — 1/2 Heisenberg model on 
(a) 4x8, (b) 8x2 and (c) 8x4 cylinders. All cylinders have 
a staggered pinning field of strength J/2 applied at the open 
ends. Because cylinder (a) is too short, edge effects push its 
bulk magnetization above the true 2D value. Cylinder (b) is 
too one-dimensional so its bulk magnetization falls below the 
2D value. Cylinder (c) has an aspect ratio of 2 : 1 which is 
almost ideal,^ making its bulk a good approximation of the 
infinite 2D system. The magnetization of cylinder (c) at the 
center rungs (m ^ 0.28) is close to the estimated 2D value of 
m 0.30743.^^ 



cylinders — the leading finite-size corrections to the en- 
ergy density of the Heisenberg model are expected to vary 
as 1/Ll7^ In Figure [e] we show the result of fitting the 
infinite cylinder energies with > 6 to this form. We 
obtain an estimate eo = —0.6694(3) consistent with the 
best published Monte Carlo result eo = -0.669437(5). ^31 

Bulk order parameters should be obtained differently 
from energies — considering infinite cylinders is not opti- 
mal. Asymptotically, infinite length cylinders of constant 
width will exhibit ID behavior. 

A particularly effective approach for measuring order 
parameters is to try to choose an optimal aspect ratio 
a = Lx/Ly for finite cylinders. The approach utilizes the 
strongly pinned boundary conditions described above. 
Say we decide to extrapolate to 2D only using cylinders 
with a fixed a. As illustrated in Figure[7) if a is too small, 
edge effects suppress fluctuations and our estimates will 
approach the 2D value from above; if a is too large the 
physics resembles that of a ID chain and our estimates 
approach from below. In several "typical" systems and 
in associated continuum theories it has been found that 
there is an optimal ratio a* where leading corrections in 
1/Ly vanish.^ Extrapolations performed at fixed a* will 
be nearly fiat, improving their accuracy. For periodic 
boundary conditions the ideal aspect ratio turns out to 
be about a* c:^ 7. For cylindrical boundaries it is much 
lower — around a* o:^ 1.9."^ This is yet another reason 
why cylindrical boundary conditions should be preferred. 
Finite size effects can be drastically reduced even if one 
does not use the exact aspect ratio. 



C. Gaps and Excited States 

In addition to ground state properties, one usually 
wants to have a good understanding of the elementary 
excitations of a system. Calculating the energy gap to 
the first excited state is of fundamental importance for 
classifying phases and estimating their robustness to per- 
turbations. The excited states themselves can be useful 
for understanding subtle orders present in the ground 
state and for computing quantities that are accessible to 
experiments. 

Traditionally there have been two ways to find excited 
states and gaps using DMRG and both methods are a 
good option for studying 2D systems. The simplest situ- 
ation is when an excited state lies in a different quantum 
number sector than the ground state. This is the case 
when calculating the spin gap of a magnet, for exam- 
ple. By taking advantage of quantum numbers DMRG 
can directly target the ground state of an excited sec- 
tor. This method is preferred as long as the excitation is 
not attracted to the open ends of the system, which can 
easily be checked from measurements of local quantities. 
If the excitation is attracted to the ends, the restricted 
sweeping approach described below is a good option. 

When studying excited states in the same sector as the 
ground state, such as singlet excitations, or when study- 
ing a model with no conserved quantities one can simulta- 
neously target the lowest few eigenstates. In this scheme 
the superblock wavefunctions for each state are kept sep- 
arate but they share one set of boundary blocks. It is 
often necessary therefore to keep an increased number 
of states in order to approximate different wavefunctions 
using a single truncated basis. 

But there are situations where neither of the above ap- 
proaches is a good option. Models of topological phases 
often have no conserved quantities plus large ground state 
degeneracy that makes multiple targeting inefficient. For 
cases like these, it is helpful to take advantage of the flexi- 
bility offered through using matrix product states (MPS) 
and matrix product operators. Representing a state as 
an MPS gives DMRG access to the entire wavefunction 
and makes it possible to compute overlaps between wave- 
functions found through separate calculations. 

This flexibility gives a new way to find excited states as 
follows. First, DMRG is used to compute a ground state 
\ipo) of the Hamiltonian H as an MPS. Then, one defines 
a Hamiltonian H' = H -\- wPq where Pq = |'0o)('^o| is a 
projection operator and w is an energy penalty for states 
not orthogonal to lipo) . w is large enough, the ground 
state iV^i) of will be the second lowest eigenstate of 
H (its first excited state or a second ground state). 

But when finding \tpi) it is neither necessary nor ad- 
visable to compute H' . One should instead work with H 
and Pq separately. During a normal Lanczos or Davidson 
step in DMRG, the product i^lV^i) gets computed in the 
current local basis. The product PolV^i) can be computed 
in a similar fashion following the procedure in Figure [S] 
The resulting tensors H\ipi) and k;Po|'0i) can then be 
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TABLE I: Comparison of Some Leading Numerical Methods for 2D Strongly Correlated Systems 



Method 


Approach 


Variational 


2D? 


Sign Problem 


Biased 


Typical Computa- 
tional Effort 


Worldline QMCP Statistical Sampling 


No 


Yes 


Yes 


No 


N 


Determinantal 


Statistical Sampling 


No 


Yes 


Yes 


No 


N^ 


Variational 
QMCP 


Statistical Sampling 


Yes 


Yes 


Fixed by Guessed 
Wavefunction 


Yes 


N^ 


Series Expansion^^^xtrapolated 
Taylor Series 


No 


Yes 


No 


Yes 


< 10 or 20 Terms 


DMRG 


Low Entanglement 


Yes 


Width <12 


No 


Very Slight^ 




PEPS 


Low Entanglement 


No^ 


Yes 


No 


Very Slight" 


D^^-D^^ (Ref.lHOl) 


MERA 


Low Entanglement 


Yes 


Yes 


No 


Very Slight"'^ 


D^^ (ReLlS]) 



a. Indicates a bias toward states having low entanglement. 

b. Though PEPS themselves are variational, observables must be computed through a controlled approximation. 

c. The pattern of tensors chosen may favor a certain type of ground state. 



of low- lying states \ipn) it is possible to obtain an even 
more accurate spectrum by computing the eigenvalues 
of H^'^ = {^l)m\H\il)n) ' This last step rotates away any 
remaining non-orthogonality present in the states \ipn) 
ensuring that energies are computed in an optimal basis. 

Each of the three methods discussed above has its mer- 
its, but one issue that affects them all is the possibility 
of spurious excitations at the open ends of the cylinder. 
In most cases such an excitation is uninteresting since it 
has no analogue in the infinite 2D system. If this prob- 
lem arises, one way to deal with it is through restricted 
sweeping. First, the ground state is computed as usual 
to very high accuracy. Then, within the same calcula- 
tion one switches over to the first excited state either by 
changing quantum numbers or through multiple target- 
ing. But now the sweeping range is restricted only to 
the bulk of the cylinder, keeping the system fixed to the 
ground state basis near the open ends. By mixing with 
this fixed basis, excitations can have tails extending be- 
yond the sweeping region, but will not get stuck in an 
unphysical state where they live only at an edge. 



IV. TWO DIMENSIONAL TENSOR NETWORK 
METHODS 

Based on DMRG's success for ID and quasi- ID sys- 
tems, it is reasonable to ask if DMRG can be extended 



added to form H'\iIji). 

Having found one can go on to compute the next 
excited state by including both Pq and Pi = 
in the effective Hamiltonian. Many low-lying states can 
be found this way with a cost that is quadratic in the 
number of states, although the method can only be 
pushed so far unless the previous wavefunctions are de- 
termined to high accuracy. Finally, after finding a set 



(a) 



(V'olV'i) 



(b) 



I 1 I 



i'olV'i) = (V'olV'i) 



V'J 



rr 



^1 



I I I 



FIG. 8: Steps for computing the product of a projector Po = 
|'0o)('0o| with a state in DMRG. After computing the 
tensors Pl and Pr that project ('0o| into the current basis, 
one (a) computes the overlap ('0o|'0i) then (b) multiplies it 
by the projected form of iV^o)- 
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to 2D in a more natural way. Early attempts to do so 
failed because they did not take into account the essential 
differences in entanglement scaling of ID and 2D systems. 
Ground states of ID gapped Hamiltonians obey an area 
law: any subsystem is entangled with its environment 
only through the boundary connecting the two. This 
means that ground states of gapped ID systems can be 
represented accurately by matrix product states with a 
fixed bond dimension. 

Most — if not all — gapped phases in 2D also obey an 
area law, but in 2D this means that the entanglement 
entropy of a subregion grows proportionally to its linear 
size. To capture this effect with an MPS requires that its 
bond dimension increase exponentially with the system 
size no matter how it is embedded into a 2D lattice. The 
key, then, to obtaining a scalable version of 2D DMRG is 
to identify new classes of variational wavefunctions that 
extend matrix product states to reproduce 2D area law 
behavior. 

Though there have been a number of promising at- 
tempts in this direction, two classes of wavefunctions 
stand out based on their usefulness in simulating a num- 
ber of realistic 2D models. The first is a natural extension 
of an MPS known as a Projected Entangled Pair State 
or PEPS.^^ Just as an MPS can be viewed as a net- 
work of rank 3 tensors, a PEPS is a network of rank 
Z + 1 tensors with Z the coordination number of the 
lattice. The usefulness of PEPS is somewhat limited by 
the steep computational cost required to optimize them, 
which grows as D^^ where D is the bond dimension. And 
although a PEPS is a well-defined variational wavefunc- 
tion, it is necessary to introduce approximations when 
computing observables, leading to energy estimates that 
are no longer variational. But because of the higher co- 
ordination number of 2D lattices, one expects that 2D 
ground states are already represented quite well by a 
PEPS with D - 10-100 versus an MPS for a ID system 
where m must be many hundreds to achieve the same 
accuracy. Perhaps the most attractive feature of the 
PEPS approach is its ability to work directly in the ther- 
modynamic limit using the so- called iPEPS method.l^^ 
The iPEPS method has been successfully applied to a 
wide variety of models including orbital models, frus- 
trated magnets, interacting bosons^^ and interacting 
fermionsl ^^ l ^-^ l 

A second type of wavefunction that has proven effective 
for 2D is the Multiscale Entanglement Renormalization 
Ansatz (MERA). Originally conceived as a concrete re- 
alization of Kadanoff 's real space renormalization group, 
a MERA generalizes an MPS by extending the ID chain 
of tensors into a layered structure where each layer rep- 
resents a coarser length scale in the RG process .1^^^ The 
MERA approach can naturally be extended to 2D by us- 
ing tensors that group patches of lattice sites together 
instead of ID segments. But unlike the PEPS ansatz, 
there is no unique MERA for a given lattice. This makes 
the method more flexible but can also be a source of 
bias if the tensors are chosen to favor a certain type of 




J \ \ \ \ \ L 

6 8 10 12 

Ly 



FIG. 9: Plots of the accuracy achievable with DMRG for 
cylinders of width Ly and a maximum number of states kept 
m. Data is for the 5* = 1/2 square lattice Heisenberg model 
with J —\. The DMRG accuracy is defined as the difference 
between the energy for a fixed m and the energy obtained 
by extrapolating in the truncation error. The horizontal lines 
show the difference between iPEPS e nerg y estimate^^ and 
the QMC estimate eo = -0.669437(5) .^Because one is al- 
ways restricted to some maximum m in DMRG, it is only 
worthwhile to extrapolate using widths Ly for which there is 
a better accuracy than iPEPS; for simulations of larger sys- 
tems it is more sensible to use PEPS or MERA directly. 



ground state. MERA optimization techniques currently 
suffer from a higher computational cost than PEPS (see 
Table |l| but a key advantage of MERA over PEPS is 
that local observables may be computed exactly and ef- 
ficiently, allowing MERA simulations to remain varia- 
tional. This property of the MERA comes from its con- 
struction in terms of isometrics that cancel outside of 
a 'light cone' emanating from the location of the oper- 
ator to be measured. Finally, like the PEPS method, 
the MERA approach can give results for both finite sys- 
tems and the thermodynamic limit. Recent applications 
of MERA to 2D systems include studies of interacting 
fermions,^^ orbital systems, topological models^ and 
frustrated magnets. 

A key question is this: when should one use 2D DMRG 
methods, and when should one try PEPS or MERA (or 
another approach)? As a first step in addressing this 
issue, we provide in Figure [9] a simple comparison be- 
tween iPEPS results'^ and DMRG data for the 2D square 
lattice Heisenberg model. (We hope that comparable 
data for MERA can be provided in the future.) The 
accuracy of iPEPS depends on the tensor dimension 
while that of DMRG depends on m and the width of 
the cylinder studied. At some width, which depends on 
how large an m or I) can be treated with current comput- 
ing resources, DMRG becomes less accurate than iPEPS. 
(Neither the values of m nor of D shown in the figure 
should be thought of as reflecting the state of the art.) 
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Whether one should use DMRG or iPEPS depends pri- 
marily on whether the bulk behavior can be extrapolated 
with cylinders smaller than this width. Alternatively, it 
can be very useful to use both DMRG and iPEPS or 
MERA, and compare results. 



We would like to acknowledge helpful discussions with 
Sasha Chernyshev, Salvatore Manmana and Philippe 
Corboz, who was kind enough to provide iPEPS data for 
the square lattice Heisenberg model. We acknowledge 
support from the NSF under DMR-0907500. 



^ Bergman D, Alicea J, Gull E, Trebst S, Balents L. 2007. 

Nature Phys. 3:487 
^ Castelnovo C, Moessner R, Sondhi S. 2008. Nature 451:42- 

45 

^ Dagotto E. 1994. Rev. Mod. Phys. 66:763-840 

^ Kastner MA, Birgeneau RJ, Shirane G, Endoh Y. 1998. 

Rev. Mod. Phys. 70:897-928 
^ Stormer HL, Tsui DC, Gossard AC. 1999. Rev. Mod. Phys. 

71:S298-S305 
^ Balents L. 2010. Nature 464:199-208 

^ Yamashita M, Nakata N, Kasahara Y, Sasaki T, Yoneyama 

N, Kobayashi N, Fujimoto S, Shibauchi T, Matsuda Y. 

2009. Nature Physics 5:44-47 
^ Yamashita M, Nakata N, Senshu Y, Nagata M, Yamamoto 

HM, Kato R, Shibauchi T, Matsuda Y. 2010. Science 

328:1246-1248 

^ Yan S, Huse DA, White SR. 2011. Science 332:1173-1176 
Giamarchi T. 2003. Quantum Physics in One Dimension. 
Oxford University Press 

Sandvik AW. 2007. Phys. Rev. Lett. 98:227202 
^2 Isakov SV, Hastings MB, Melko RG. arxiv: 1102. 1721 
White SR, Scalapino DJ, Sugar RL, Loh EY, Gubernatis 
JE, Scalettar RT. 1989. Phys. Rev. B 40:506-516 
Troyer M, Wiese UJ. 2005. Phys. Rev. Lett. 94:170201 
White SR. 1992. Phys. Rev. Lett. 69:2863-2866 
White S. 1993. Phys. Rev. B 48:10345-10356 
^'^ Liischer A, LauchU AM. 2009. Phys. Rev. B 79:195102 
Poilblanc D, Mambrini M, Schwandt D. 2010. Phys. Rev. 
B 81:180402(R) 

White SR, Chernyshev A. 2007. Phys. Rev. Lett. 99:127004 
2° Schollwock U. 2005. Rev. Mod. Phys. 77:259-315 
2^ Hastings MB. 2007. J. Stat. Mech. 2007:P08024 
2^ White S, Huse D. 1993. Phys. Rev. B 48:3844-3852 
2^ Schollwock U, White SR. cond-mat/ 06060 18 
2^ White SR. 2009. Phys. Rev. Lett. 102:190601 

Stoudenmire E, White SR. 2010. New Journal of Physics 

12:32 

2^ Ostlund S, Rommer S. 1995. Phys. Rev. Lett. 75:3537-3540 
^'^ Jiang HC, Weng ZY, Xiang T. 2008. Phys. Rev. Lett. 
101:090603 

2^ Evenbly G, Vidal G. 2010. Phys. Rev. Lett. 104:187203 

Wang L, Pizorn I, Verstraete F. 2011. Phys. Rev. B 

83:134421 
^° McCulloch I. arxiv:0804.25090 

Kallin AB, Hastings MB, Melko RG, Singh RRP. 

arxiv:1107.2840 

Schollwock U. 2011. Annals of Physics 326:96-192 
Hallberg KA. 2006. Advances in Physics 55:477-526 
Noack RM, Manmana SR. 2005. AIP Conference Proceed- 
ings 789:93-163 

White SR. 1996. Phys. Rev. Lett. 77:3633-3636 
Yoshikawa S, Okunishi K, Senda M, Miyashita S. 2004. J. 
Phys. Soc. Jpn. 73:1798-1804 

Weng M, Sheng DN, Weng ZY, Bursill RJ. 2006. Phys. 



Rev. B 74 

Jiang HC, Weng ZY, Sheng DN. 2008. Phys. Rev. Lett. 
101:117203 

Capriotti L, Scalapino DJ, White SR. 2004. Phys. Rev. B 
93:177004 

^° Jiang HC, Kriiger F, Moore JE, Sheng DN, Zaanen J, 
Weng ZY. 2009. Phys. Rev. B 79:174409 
Sheng DN, Motrunich 01, Fisher MPA. 2009. Phys. Rev. 
B 79:205112 

^2 Block MS, Sheng DN, Motrunich 01, Fisher MPA. 2011. 
Phys. Rev. Lett. 106:157202 

Jiang HC, Gu ZC, Qi XL, Trebst S. 2011. Phys. Rev. B 
83:245104 

White SR, Noack RM, Scalapino DJ. 1994. Phys. Rev. Lett. 

73:886-889 

du Croo de Jongh M, van Leeuwen J. 1998. Phys. Rev. B 
57:8494-8500 

Zhao J, Wang X, Xiang T, Su Z, Yu L, Lou J, Chen C. 
2006. Phys. Rev. B 73:012411 

White SR, Scalapino DJ. 1997. Phys. Rev. B 55:R14701- 
R14704 

White SR, Scalapino DJ. 1998. Phys. Rev. Lett. 80:1272- 
1275 

White SR, Scalapino DJ. 1998. Phys. Rev. Lett. 81:3227- 
3230 

^° Tohyama T, Gazza C, Shih CT, Chen YC, Lee TK, 
Maekawa S, Dagotto E. 1999. Phys. Rev. B 59:R11649- 
R11652 

White SR, Scalapino DJ. 1999. Phys. Rev. B 60:R753- 
R756 

^2 Kampf AP, Scalapino DJ, White SR. 2001. Phys. Rev. B 

64:052509 

Emery VJ, Kivelson SA, Lin HQ. 1990. Phys. Rev. Lett. 

64:475-478 

White SR, Scalapino DJ. 2000. Phys. Rev. B 61:6320-6326 
Scalapino DJ, White SR. 2001. Foundations of Physics 
31:27-39 

Chernyshev AL, White SR, Neto AHC. 2002. Physica B: 
Condensed Matter 312-313:566 - 568 

Chernyshev AL, White SR, Castro Neto AH. 2002. Phys. 
Rev. B 65:214527 

White SR, Scalapino DJ. 2004. Phys. Rev. B 70:220506(R) 

Chernyshev AL, Castro Neto AH, White SR. 2005. Phys. 

Rev. Lett. 94:036407 
^° White SR, Scalapino DJ. 2009. Phys. Rev. B 79:220504 

White SR, Scalapino DJ. 2003. Phys. Rev. Lett. 91:136403 
^2 Hager G, Wellein G, Jeckelmann E, Fehske H. 2005. Phys. 

Rev. B 71:075108 

Fehske H, Hager G, Wellein G, Jeckelmann E. 2006. Phys- 
ica B: Cond. Mat. 378-380:319 - 320 

Jiang HC, Weng M, Weng ZY, Sheng DN, Balents L. 2009. 
Phys. Rev. B 79 

Motrunich 01, Fisher MPA. 2007. Phys. Rev. B 75:235116 
Sheng DN, Motrunich 01, Trebst S, Gull E, Fisher MPA. 



11 



2008. Phys. Rev. B 78:054520 

Feiguin AE, Fisher MPA. 2011. Phys. Rev. B 83:115104 
Nishimoto S, Nakamura M, O'Brien A, Fulde P. 2010. 
Phys. Rev. Lett. 104:196401 

Jeckelmann E, White SR. 1998. Phys. Rev. B 57:6376-6385 
^° Liang S, Pang H. 1994. Phys. Rev. B 49:9214-9217 

Dukelsky J, Martin-Delgado MA, Nishino T, Sierra G. 
1998. Europhys. Lett. 43:457 
^2 White SR. 2005. Phys. Rev. B 72:180403 

Neuberger H, Ziman T. 1989. Phys. Rev. B 39:2608-2618 
Sandvik AW. 1997. Phys. Rev. B 56:11678-11690 
Sandvik AW, Evertz HG. 2010. Phys. Rev. B 82:024407 
Sandvik AW. 2010. AIP Conference Proceedings 1297:135- 
338 

Blankenbecler R, Scalapino DJ, Sugar RL. 1981. Phys. 
Rev. D 24:2278-2286 
'^^ Foulkes WMC, Mitas L, Needs RJ, Rajagopal G. 2001. 
Rev. Mod. Phys. 73:33-83 

Domb C, Green M. 1974. Phase Transitions and Critical 
Phenomena, vol. 3. Academic Press, London 
Corboz P. (private communication) 

Corboz P, Evenbly G, Verstraete F, Vidal G. 2010. Phys. 
Rev. A 81:010303 



Nishino T, Hieida Y, Okunishi K, Maeshima N, Akutsu Y, 
Gendiar A. 2001. Prog. Theor. Phys. 105:409-417 
Verstraete F, Porras D, Cirac JL 2004. Phys. Rev. Lett. 
93:227205 

Verstraete F, Murg V, Cirac JL 2008. Advances in Physics 
57:143-224 

Jordan J, Orus R, Vidal G, Verstraete F, Cirac JL 2008. 

Phys. Rev. Lett. 101:250602 

Orus R, Vidal G. 2009. Phys. Rev. B 80:094403 

Orus R, Doherty AC, Vidal G. 2009. Phys. Rev. Lett. 

102:077203 

Bauer B, Vidal G, Troyer M. 2009. J. Stat. Mech 
2009:P09006 

Jordan J, Oriis R, Vidal G. 2009. Phys. Rev. B 79:174515 
^° Corboz P, Orus R, Bauer B, Vidal G. 2010. Phys. Rev. B 
81:165104 

Corboz P, Jordan J, Vidal G. 2010. Phys. Rev. B 82:245119 
^2 Vidal G. 2007. Phys. Rev. Lett. 99:220405 

Evenbly G, Vidal G. 2009. Phys. Rev. B 79:144108 
Cincio L, Dziarmaga J, Oles AM. 2010. Phys. Rev. B 
82:104416 

Tagliacozzo L, Vidal G. 2011. Phys. Rev. B 83:115127 



